Benthic Dives Investigations

Authors

Astarte Brown

Conner Hale

Joffrey Joumaa

Published

March 29, 2023

1 Import library

Code
# data viz
library(ggplot2)
library(ggh4x)
library(ggOceanMaps)
library(patchwork)
library(viridis)
library(ggpubr)
library(grid)
library(viridis)
library(ggnewscale)

# table
library(gt)

# map
library(rnaturalearth)
library(ggmap)
library(ggsn)
library(sf)
library(sp)
library(smoothr)

# kernel density
library(eks)
library(ggdensity)

# stat
library(Hmisc)
library(circular)
library(CircStats)

# char
library(stringr)

# solar angle
library(oce)

# data wrangling
library(tidyr)
library(dplyr)
library(purrr)
library(forcats)
library(lubridate)

# model
library(lme4)
library(fitdistrplus)

2 Setting up custom function

2.1 windrose

Code
windrose <-
  function(data_to_plot,
           grid = NULL,
           set_title = NULL,
           legend_position = "none",
           facet = F) {
    # this code comes from Roxanne
    uniqhours <- 1:24 * (360 / 24)

    # trick to align hours om the graph
    data_to_plot <- rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])

    for (i in 1:24) {
      # turn hours to radians
      if (i == 1) {
        temp <- rep(uniqhours[i], data_to_plot$nb_ind_hour[i])
        day_night <- rep("night", data_to_plot$nb_ind_hour[i])
      } else {
        temp <- c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i]))
        day_night <- c(day_night, rep(
          if_else(between(i, 7, 20), "day", "night"),
          data_to_plot$nb_ind_hour[i]
        ))
      }
    }
    data2 <- data.frame(direction = temp)

    deg <- 15 # choose bin size (degrees/bin)
    dir.breaks <- seq(0 - (deg / 2), 360 + (deg / 2), deg) # define the range of each bin

    # assign each direction to a bin range
    dir.binned <-
      cut(data2$direction,
        breaks = dir.breaks,
        ordered_result = TRUE
      )
    # generate pretty lables
    dir.labels <- as.character(c(seq(0, 360 - deg, by = deg), 0))

    # replace ranges with pretty bin lables
    levels(dir.binned) <- dir.labels

    # Assign bin names to the original data set
    data2$dir.binned <- dir.binned

    # add origin if any
    if (facet) {
      data2$origin <- unique(data_to_plot$origin)
    }

    # set up max value
    maxvalue <- 35

    # initialise the plot
    plt.dirrose_2 <- ggplot()

    # check if grid
    if (!is.null(grid)) {
      plt.dirrose_2 <- plt.dirrose_2 +
        geom_hline(
          yintercept = grid,
          colour = "grey20",
          linewidth = .2
        )
    }
    plt.dirrose_2 <- plt.dirrose_2 +
      geom_vline(
        xintercept = c(seq(1, 24, 2)),
        colour = "grey30",
        linewidth = 0.2
      ) + # 24 vertical lines at center of the 30? ranges.
      geom_hline(
        yintercept = maxvalue,
        colour = "white",
        linewidth = .5
      ) + # Darker horizontal line as the top border (max).
      # On top of everything we place the histogram bars.
      geom_bar(
        data = data2,
        aes(x = dir.binned, fill = day_night),
        width = 1,
        colour = "black",
        linewidth = 0.3
      ) +
      # geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +
      scale_x_discrete(
        drop = FALSE,
        labels = c(
          0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, ""
        )
      ) +
      scale_fill_manual(values = c("white", "darkgrey"), name = "Time of day") +
      labs(x = "Time (hours)", y = "Count", title = set_title) +
      coord_polar(start = -(deg / 2) * (pi / 180))

    # if facet
    if (facet) {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(. ~ origin)
    }

    # Wraps the histogram into a windrose
    plt.dirrose_2 <- plt.dirrose_2 +
      theme_bw() +
      theme(
        legend.position = legend_position,
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key.size = unit(0.5, "cm")
      )

    # return
    return(plt.dirrose_2)
  }

2.2 getNightLength

Based on here: https://stackoverflow.com/questions/59566553/calculate-night-length-in-r

Code
getNightLength <- function(survey_date, latitude, longitude) {
  # survey_date <- strptime(survey_date, "%d/%m/%y")
  time_origin <- strptime("01/01/2000 12:00", "%d/%m/%Y %H:%M")
  n <- as.numeric(difftime(survey_date, time_origin, units = "days"))
  Jstar <- n - longitude / 360
  M <- (357.5291 + 0.98560028 * Jstar) %% 360
  C <- 1.9148 * sin(M * 2 * pi / 360) +
    0.0200 * sin(2 * M * 2 * pi / 360) +
    0.0003 * sin(3 * M * 2 * pi / 360)
  lambda <- (M + C + 180 + 102.9372) %% 360
  Jtransit <- as.double(2451545.000) +
    as.double(Jstar) +
    as.double(0.0053 * sin(M * 2 * pi / 360)) -
    as.double(0.0069 * sin(lambda * 4 * pi / 360))
  sindelta <- sin(lambda * 2 * pi / 360) * sin(23.44 * 2 * pi / 360)
  delta <- asin(sindelta) * 360 / (2 * pi)
  cos_omega <- (sin(-0.83 * 2 * pi / 360) - (sindelta * sin(latitude * 2 * pi / 360))) /
    (cos(latitude * 2 * pi / 360) * cos(delta * 2 * pi / 360))
  omega <- acos(cos_omega) * 360 / (2 * pi)
  Jrise <- Jtransit - omega / 360
  Jset <- Jtransit + omega / 360
  return((1 - (Jset - Jrise)) * 24)
}

3 Import data

Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.

Code
# import the data
data_dive <- readRDS("../export/data_dive.rds")

# filter on seals departing from ANNU
data_dive <- data_dive %>%
  filter(DepartureLocation == "ANNU")
Code
# get solar angle when data location
solar_angle_inter <- data_dive %>%
  filter(!is.na(Lat)) %>%
  mutate(solar_angle = sunAngle(date, Long, Lat)$altitude) %>%
  select(id, DiveNumber, solar_angle)

# merge that with data_dive
data_dive <- merge(data_dive,
  solar_angle_inter,
  by = c("id", "DiveNumber"),
  all.x = T
)

# add day-night
# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/
data_dive <- data_dive %>%
  mutate(day_night = if_else(solar_angle < -18, "night", "day"))

Let’s add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, i.e. if this difference is above xxx m then is not considered as a benthic dive.

Code
data_dive <- data_dive %>%
  # change DiveType to include condition on difference
  # between max depth and bathymetry (<50m)
  mutate(
    DiveType_50 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 50, 0, DiveType),
    DiveTypeName_50 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 50,
      "Transit",
      DiveTypeName
    ),
    # same with <100m
    DiveType_100 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 100, 0, DiveType),
    DiveTypeName_100 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 100,
      "Transit",
      DiveTypeName
    ),
    # same with <150m
    DiveType_150 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 150, 0, DiveType),
    DiveTypeName_150 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 150,
      "Transit",
      DiveTypeName
    ),
    # same with <200m
    DiveType_200 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 200, 0, DiveType),
    DiveTypeName_200 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 200,
      "Transit",
      DiveTypeName
    )
  )

4 Maps

Code
# check if background_ggoceanmaps exist
if (file.exists("../export/background_ggoceanmap.rds")) {
  trip <- readRDS("../export/background_ggoceanmap.rds")
} else {
  # using ggOceanMaps
  trip <- basemap(
    limits = c(170, -110, 20, 59),
    bathymetry = TRUE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip <- trip +
    geom_sf(data = graticule, color = "grey50")

  trip$layers <- trip$layers[c(1, 3, 2)]

  # save result
  saveRDS(trip, "../export/background_ggoceanmap.rds")
}
Code
# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(
  x = df_kernel_dens_sf,
  # boundary = ne_coastline(scale = 110, returnclass = "sf"),
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  # ggplot(data=ne_coastline(scale = 110, returnclass = "sf"))+geom_sf()+
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # kde_test,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "All Benthic Dives",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 1: Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code
# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_50 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<50",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 2: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code
# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_100 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<100",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 3: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code
# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_150 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<150",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 4: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code
# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_200 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<200",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )


# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 5: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

5 Windrose

Code
# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(3000, 6000, 9000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 6: Circular histogram plots displaying the time when they perform benthic dives across their foraging trip.

Code
# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.5 (p-value: 0)"
Code
# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 7: Circular histogram plots displaying the time when they perform benthic dives (diff<50) across their foraging trip.

Code
# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_50 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
Code
# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 8: Circular histogram plots displaying the time when they perform benthic dives (diff<100) across their foraging trip.

Code
# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_100 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
Code
# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 9: Circular histogram plots displaying the time when they perform benthic dives (diff<150) across their foraging trip.

Code
# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_150 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
Code
# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 10: Circular histogram plots displaying the time when they perform benthic dives (diff<200) across their foraging trip.

Code
# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_200 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"

6 Boxplot of the benthic dive rate (day vs. night)

6.0.1 Graph

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 11: Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time

6.0.2 Model

Code
# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)
summary statistics
------
min:  1   max:  4.333333 
median:  1.555556 
mean:  1.669479 
estimated sd:  0.4221527 
estimated skewness:  1.473928 
estimated kurtosis:  6.414337 

Figure 12: Result from descdist function to determine the best distribution to use for data modelling

Based on Figure 12, it seems that a gamma distribution would be the best to model this dataset.

Code
# model
mdl_all <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_all) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 0.64660385 0.005864671 110.254069 0.000000e+00
fixed NA day_nightday -0.02527306 0.004559170 -5.543347 2.967447e-08
fixed NA scale(mean_dive_duration) 0.11434250 0.004412862 25.911191 4.981557e-148
ran_pars id sd__(Intercept) 0.06307474 NA NA NA
ran_pars Residual sd__Observation 0.11146813 NA NA NA

6.0.3 Go Further…

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

Code
# model
mdl_all_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_all_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 5.41274751 0.041968504 128.97166 0.000000e+00
fixed NA day_nightnight -0.72679883 0.006652905 -109.24533 0.000000e+00
fixed NA scale(mean_dive_duration) 0.27370713 0.007918629 34.56497 8.494370e-262
fixed NA scale(mean_period_duration) 0.08335783 0.002524341 33.02162 3.975539e-239
ran_pars id sd__(Intercept) 0.83810984 NA NA NA
Code
qqnorm(residuals(mdl_all_extra))
qqline(residuals(mdl_all_extra))

6.0.4 Graph

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 13: Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time

6.0.5 Model

Code
# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)
summary statistics
------
min:  1   max:  5.166667 
median:  2.555556 
mean:  2.581992 
estimated sd:  0.6479778 
estimated skewness:  0.3352311 
estimated kurtosis:  3.477161 

Figure 14: Result from descdist function to determine the best distribution to use for data modelling

Based on Figure 14, it seems that a gamma distribution would be the best to model this dataset.

Code
# model
mdl_50 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_50) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 0.41544372 0.003997844 103.916932 0.000000e+00
fixed NA day_nightday -0.01594999 0.003819268 -4.176192 2.964297e-05
fixed NA scale(mean_dive_duration) 0.08109884 0.002975447 27.256018 1.410142e-163
ran_pars id sd__(Intercept) 0.03423312 NA NA NA
ran_pars Residual sd__Observation 0.14226941 NA NA NA

6.0.6 Go Further…

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

Code
# model
mdl_50_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_50_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 4.0982271 0.042556997 96.29972 0.000000e+00
fixed NA day_nightnight -0.5263318 0.008752889 -60.13235 0.000000e+00
fixed NA scale(mean_dive_duration) -0.1472508 0.011465763 -12.84265 9.457963e-38
fixed NA scale(mean_period_duration) 0.1262806 0.004555875 27.71820 4.213980e-169
ran_pars id sd__(Intercept) 0.8416426 NA NA NA
Code
qqnorm(residuals(mdl_50_extra))
qqline(residuals(mdl_50_extra))

6.0.7 Graph

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 15: Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time

6.0.8 Model

Code
# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)
summary statistics
------
min:  1   max:  5 
median:  2.565217 
mean:  2.612124 
estimated sd:  0.643816 
estimated skewness:  0.4024764 
estimated kurtosis:  3.504268 

Figure 16: Result from descdist function to determine the best distribution to use for data modelling

Based on Figure 16, it seems that a gamma distribution would be the best to model this dataset.

Code
# model
mdl_100 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_100) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 0.41064106 0.003745655 109.631316 0.000000e+00
fixed NA day_nightday -0.01556972 0.003443919 -4.520931 6.156812e-06
fixed NA scale(mean_dive_duration) 0.07981819 0.002733776 29.197048 2.113987e-187
ran_pars id sd__(Intercept) 0.03328018 NA NA NA
ran_pars Residual sd__Observation 0.13201902 NA NA NA

6.0.9 Go Further…

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

Code
# model
mdl_100_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_50_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 4.0982271 0.042556997 96.29972 0.000000e+00
fixed NA day_nightnight -0.5263318 0.008752889 -60.13235 0.000000e+00
fixed NA scale(mean_dive_duration) -0.1472508 0.011465763 -12.84265 9.457963e-38
fixed NA scale(mean_period_duration) 0.1262806 0.004555875 27.71820 4.213980e-169
ran_pars id sd__(Intercept) 0.8416426 NA NA NA
Code
qqnorm(residuals(mdl_100_extra))
qqline(residuals(mdl_100_extra))

6.0.10 Graph

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 17: Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time

6.0.11 Model

Code
# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)
summary statistics
------
min:  1   max:  5 
median:  2.565217 
mean:  2.610992 
estimated sd:  0.6399849 
estimated skewness:  0.4233921 
estimated kurtosis:  3.460216 

Figure 18: Result from descdist function to determine the best distribution to use for data modelling

Based on Figure 18, it seems that a gamma distribution would be the best to model this dataset.

Code
# model
mdl_150 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_150) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 0.41135058 0.003639353 113.0285 0.000000e+00
fixed NA day_nightday -0.01641792 0.003291418 -4.9881 6.097591e-07
fixed NA scale(mean_dive_duration) 0.08012749 0.002635001 30.4089 4.188961e-203
ran_pars id sd__(Intercept) 0.03289909 NA NA NA
ran_pars Residual sd__Observation 0.12728439 NA NA NA

6.0.12 Go Further…

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

Code
# model
mdl_150_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_150_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 4.1866975 0.042103917 99.43725 0.000000e+00
fixed NA day_nightnight -0.5254813 0.008305595 -63.26835 0.000000e+00
fixed NA scale(mean_dive_duration) -0.1149068 0.010908400 -10.53379 6.035410e-26
fixed NA scale(mean_period_duration) 0.1365982 0.004237414 32.23623 5.487327e-228
ran_pars id sd__(Intercept) 0.8334743 NA NA NA
Code
qqnorm(residuals(mdl_150_extra))
qqline(residuals(mdl_150_extra))

6.0.13 Graph

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 19: Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time

6.0.14 Model

Code
# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)
summary statistics
------
min:  1   max:  5 
median:  2.541667 
mean:  2.600331 
estimated sd:  0.6434467 
estimated skewness:  0.3851526 
estimated kurtosis:  3.429033 

Figure 20: Result from descdist function to determine the best distribution to use for data modelling

Based on Figure 20, it seems that a gamma distribution would be the best to model this dataset.

Code
# model
mdl_200 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_200) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 0.41345050 0.003696624 111.845416 0.000000e+00
fixed NA day_nightday -0.01702406 0.003423988 -4.971998 6.626632e-07
fixed NA scale(mean_dive_duration) 0.08153323 0.002723215 29.940067 5.926044e-197
ran_pars id sd__(Intercept) 0.03261119 NA NA NA
ran_pars Residual sd__Observation 0.13047509 NA NA NA

6.0.15 Go Further…

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

Code
# model
mdl_200_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_200_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")
GLMM
effect group term estimate std.error statistic p.value
fixed NA (Intercept) 4.2052051 0.042347958 99.30125 0.000000e+00
fixed NA day_nightnight -0.5284598 0.008190191 -64.52351 0.000000e+00
fixed NA scale(mean_dive_duration) -0.1162468 0.010784403 -10.77916 4.318092e-27
fixed NA scale(mean_period_duration) 0.1405915 0.004144919 33.91899 3.496541e-252
ran_pars id sd__(Intercept) 0.8386222 NA NA NA
Code
qqnorm(residuals(mdl_200_extra))
qqline(residuals(mdl_200_extra))

7 Boxplot of the number of benthic dive (day vs. night)

Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 1000)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 21: Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time

Code
broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")
nb_benthic_dive ~ day_night
term estimate std.error statistic p.value
(Intercept) 5.96 0.05 117.16 0.00
day_nightnight −0.84 0.09 −9.07 0.00
Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 22: Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time

Code
broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")
nb_benthic_dive ~ day_night
term estimate std.error statistic p.value
(Intercept) 4.54 0.07 65.08 0.00
day_nightnight −0.54 0.12 −4.69 0.00
Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 23: Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time

Code
broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")
nb_benthic_dive ~ day_night
term estimate std.error statistic p.value
(Intercept) 4.61 0.07 65.55 0.00
day_nightnight −0.54 0.12 −4.65 0.00
Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 24: Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time

Code
broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")
nb_benthic_dive ~ day_night
term estimate std.error statistic p.value
(Intercept) 4.64 0.07 65.52 0.00
day_nightnight −0.54 0.12 −4.63 0.00
Code
data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 25: Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time

Code
broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")
nb_benthic_dive ~ day_night
term estimate std.error statistic p.value
(Intercept) 4.67 0.07 65.47 0.00
day_nightnight −0.55 0.12 −4.66 0.00

8 Extra

Let’s look at the distribution of the dive duration to decide a threshold.

Code
data_dive %>%
  ggplot(aes(x = Dduration)) +
  geom_histogram(
    fill = "grey80",
    color = "black"
  ) +
  geom_vline(xintercept = 23 * 60, linetype = "dashed", col = "lightgreen") +
  geom_vline(xintercept = 4000, linetype = "dashed", col = "salmon") +
  labs(y = "# of dives", x = "Dive duration (s)") +
  theme_minimal() +
  theme(axis.line = element_line(arrow = arrow(type = "closed", length = unit(0.1, "inches"))))

Ok, so let’s pick 2500 sec and look on a map where this dives are.

Code
# define palette
data_inter <- data_dive %>%
  filter(!is.na(Lat) &
    Dduration > 4000) %>%
  arrange(id, date) %>%
  mutate(
    date_post = date + Dduration,
    diff_time = as.numeric(
      c(difftime(
        tail(date, -1),
        head(date_post, -1)
      ), 0),
      units = "mins"
    )
  )

colPal <- scales::col_numeric(
  palette = "RdYlGn",
  domain = c(0, max(data_inter$diff_time))
)
trip +
  # new scale
  new_scale_fill() +
  # kernel
  ggspatial::geom_spatial_point(
    data = data_inter,
    aes(x = Long, y = Lat),
    crs = 4326,
    col = colPal(data_inter$diff_time),
    size = 1
  ) +
  # no display alpha
  guides() +
  # title
  labs(
    title = "Dives longer than 4000 seconds",
    subtitle = paste0(nrow(data_inter), " dives")
  )